NASA Contractor Report 4519 


e- 7 ^/ 


An Analysis for High Reynolds Number 
Inviscid/Viscid Interactions in Cascades 


Mark Barnett, Joseph M. Verdon, and Timothy C. Ayer 
United Technologies Research Center 
East Hartford, Connecticut 


May 1993 


Prepared for the 
Lewis Research Center 
Under Contract NAS3-25425 


IW\S/\ 


An Analysis for High Reynolds Number 
Inviscid/ Viscid Interactions in Cascades 


Contents 

Summary 1 

1 Introduction 2 

2 General Concepts 4 

3 Inviscid Region 6 

4 Viscous Layer 8 

5 Inviscid/Viscid Interaction 12 

6 Numerical Examples 13 

6.1 Compressor Exit Guide Vane 14 

6.2 High-Speed Compressor 15 

6.3 Turbine 16 

6.4 Timing Study and Convergence Behavior 17 

7 Concluding Remarks 19 

References 20 

List of Symbols 23 

List of Figures 26 

Figures 1 through 13 27 

A Details of the Viscous-Layer Solution Procedure 40 

A.l Finite-Difference Approximations 40 

A. 2 Surface Procedure 42 

A. 3 Wake Procedure 44 


l 


AN ANALYSIS FOR HIGH REYNOLDS NUMBER INVISCID/VISCID 


INTERACTIONS IN CASCADES* 


Mark Barnett, Joseph M. Verdon, and Timothy C. Ayer 
United Technologies Research Center 
East Hardford, Connecticut 


Summary 


An efficient steady analysis for predicting strong inviscid/viscid interaction phenomena 
such as viscous-layer separation, shock/boundary-layer interaction and trailing-edge/near- 
wake interaction in turbomachinery blade passages is needed as part of a comprehensive 
analytical blade design prediction system. Such an analysis is described in the present re- 
port. It uses an inviscid/viscid interaction approach, in which the flow in the outer inviscid 
region is assumed to be potential, and that in the inner or viscous-layer region is governed 
by Prandtl’s equations. The inviscid solution is determined using an implicit, least-squares, 
finite-difference approximation, the viscous-layer solution using an inverse, finite-difference, 
space-marching method which is applied along the blade surfaces and wake streamlines. 
The inviscid and viscid solutions are coupled using a semi-inverse global iteration procedure, 
which permits the prediction of boundary-layer separation and other strong-interaction phe- 
nomena. Results are presented for three cascades, with a range of inlet flow conditions 
considered for one of them, including conditions leading to large-scale flow separations. 
Comparisons with Navier-Stokes solutions and experimental data are also given. 


*Prepared for Lewis Research Center under Contract NAS3-25425. 


1. Introduction 


An important problem faced by engine designers is the prediction of high Reynolds num- 
ber (Re) viscous flow and, in particular, viscous separation phenomena in compressor and 
turbine blade passages. Viscous effects control aerodynamic losses, heat transfer rates and 
stall, and hence, must be accounted for. In addition to steady-flow applications, the ability 
to account for viscous effects in unsteady flows is needed for aeroelastic and aeroacoustic 
design predictions, e.g., to predict the onset of stall flutter, blade row interactions due to the 
convection of viscous wakes from upstream rows and other unsteady effects which impact 
the structural and acoustic characteristics of turbomachinery blade rows. Clearly, efficient 
analytical procedures for predicting steady and unsteady viscous flows in high-performance 
compressor and turbine blading would be a significant contribution to a successful blade 
design prediction system. 

The analysis described in this report is being developed as part of a research program 
which has the goal of constructing reliable and efficient theoretical prediction methods for 
steady and unsteady viscous flows, at high Reynolds numbers, through two- and quasi- 
three-dimensional subsonic and transonic cascades. The approach to be followed is similar 
to that which has been applied successfully in external aerodynamics, where inviscid/viscid 
interaction (IVI) concepts have been used to predict the complete steady (e.g., Refs. [1] and 
[2]), and unsteady (e.g., Refs. [3] and [4]) flow fields. 

Thus, for example, for high Reynolds number flow through a cascade, the complete flow 
can be divided conceptually into two regions: an “outer” inviscid region and an “inner” 
viscous region. The construction of a general viscous cascade solver involves first, the de- 
velopment of component flow solvers, and second, the implementation of these component 
solvers into an overall computational procedure to produce a complete analysis. Solution 
methods for steady subsonic and transonic inviscid flows through cascades (e.g., Ref. [5]) 
and for steady boundary-layer and wake flows (e.g., Refs. [6] and [7]) have been developed 
to a relatively mature state. Methods for coupling such solutions have also been developed 
and assessed through a number of model problem studies (e.g., Refs. [6]-[8]). Inviscid/viscid 
interaction procedures for predicting steady flow in cascades have been developed (e.g., 
Refs. [9]— [13] ) and applied over a wide range of inlet flow conditions, including conditions 
leading to stall [13]. A similar approach can be developed for unsteady flows; for example, 
a linearized inviscid analysis [14] and a weak- interact ion unsteady viscous-layer analysis [15] 
have been coupled and applied to predict attached flows [16]. 

The focus of this report is on the development of an accurate and efficient steady cascade 
analysis that will provide the foundation for an unsteady procedure to be developed later. 
The ability of a steady inviscid/viscid interaction analysis to predict compressor cascade 
performance with accuracy comparable to that of a state-of-the-art Navier-Stokes (N-S) 
analysis has been demonstrated in Ref. [13]. Here, an Euler inviscid calculation was coupled 
to a finite-difference viscous-layer analysis. Although accurate, this approach, like current 
N-S analyses, requires too much computational effort for repetitive design calculations. The 
present effort is intended to realize the additional goal, beyond that of accuracy, of providing 
an efficient analysis suitable for design applications. 

In the present approach, the “outer” inviscid flow is determined using potential-flow the- 
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ory; thus, it is assumed that the inviscid flow is isentropic and irrotational. The steady, 
cascade, full-potential analysis (SFLOW) of Hoyniak and Verdon [17] is employed. SFLOW 
was constructed for use with the LINearized inviscid unsteady FLOw (LINFLO) analysis of 
Verdon and Caspar [14], to provide a comprehensive and compatible steady and unsteady in- 
viscid flow prediction capability for cascades. In the present calculation procedure (which will 
be referred to as SFLOW-IVI), viscous effects are incorporated by adjusting the blade and 
wake surface boundary conditions in SFLOW to account for the effects of viscous displace- 
ment thickness. The nonlinear inviscid analysis, coupled with the IVI iteration procedure, 
allows nonlinear changes to the base flow due to viscous effects to be evaluated. The ability 
to treat nonlinear perturbations is especially important in transonic flows in which shock 
positions are significantly altered by viscous displacement effects. Although the analysis 
described in this report is presently restricted to subsonic flows, it will be extended to treat 
transonic flows in the future. 
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2. General Concepts 


For flows of practical interest in either external or internal aerodynamics, the Reynolds 
number is usually sufficiently high so that the flow past an airfoil or blade can be divided 
into two regions: an “inner” dissipative region consisting of boundary layers and wakes, and 
an “outer” inviscid region. The principal interaction between the flows in the viscid and 
inviscid regions arises from the displacement thickness effect which leads to thickened semi- 
infinite equivalent bodies with corresponding changes in surface pressures. If the interaction 
is “weak,” then the complete flow problem can be solved sequentially. This traditional 
approach for calculating the interaction between the inviscid and viscous parts of the flow 
is based on a direct hierarchy between the two regions, which is applicable as long as the 
disturbances to the inviscid flow due to the viscous displacement effect remain weak. 

Flows over airfoils, however, involve both a weak overall interaction arising from standard 
displacement thickness and wake curvature effects, and local strong-displacement interac- 
tions caused, for example, by viscous-layer separations, shock/boundary-layer interactions, 
and trailing-edge/near-wake interactions. These features can lead to singularities in a clas- 
sical boundary- layer solution and a subsequent breakdown of a weak-interaction solution 
procedure. In addition, viscous displacements in the strong-interaction region cause sub- 
stantial changes in the local inviscid pressure field and can, in some cases (e.g., in flows with 
large-scale separations), cause substantial changes in the global pressure field as well. The 
concept of an inner viscous region and an outer inviscid region still holds, but the classi- 
cal hierarchical structure of the flow no longer applies. Thus, in a local strong-interaction 
region, the hierarchy changes from “direct” (i.e., pressure determined by the inviscid flow) 
to “interactive” (i.e., pressure determined by a mutual interaction between the inviscid and 
the viscous-layer flows), and this change must be accommodated within a comprehensive 
inviscid/ viscid interaction analysis. 

In a classical weak-interaction calculation, the blade surface and wake streamline pres- 
sure distribution is determined from the zeroth-order inviscid solution, that is, the inviscid 
solution unperturbed by viscous effects. This distribution is imposed when solving the 
viscous-layer equations using the direct approach mentioned above. The resulting displace- 
ment thickness distribution, 6 ( 5 ), is then used to obtain the first-order (perturbation) inviscid 
solution, thereby accounting for the changes to the base inviscid flow due to viscous displace- 
ment effects. The resulting changes in the blade pressure distributions and in the downstream 
freestream flow properties (e.g., Mach number and flow angle) can then be calculated. It is 
sometimes possible to continue this sequential solution procedure until a converged solution 
for the entire flow is achieved. This hierarchical approach, known as Prandtl iteration, works 
well if the interaction between the flows in the viscid and inviscid regions is weak. However, 
local regions of strong inviscid/viscid interaction are inevitably present in realistic airfoil 
flows. When the strong interaction leads to viscous-layer separation the Prandtl iteration 
usually fails and an alternative approach must be used. 

It is well known that the difficulties associated with the Prandtl iteration stem from 
the specification of the pressure in the viscous layer. If the boundary layer separates, the 
Goldstein singularity [18] prevents a continuation of the direct viscous-layer solution down- 
stream of the predicted separation point [19]. This difficulty can be circumvented by using 
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an inverse viscous-layer solution procedure in which the displacement thickness is specified 
instead of the pressure. The inverse procedure permits viscous solutions to be continued 
through local strong-interaction regions, including separated flow regions. 

The approach taken here employs an I VI model to calculate high Reynolds number [Re) 
flows through two-dimensional cascades. The flow in the outer inviscid region is governed 
by the full-potential equation and that in the inner viscous region by Prandtl’s viscous- 
layer equations. The non-hierarchical nature of strong interactions is accounted for in the 
procedure used to couple the two solutions. 

We consider high Reynolds number [Re = p* V'^L* / /i*^) steady flow, with negligible 
body forces, of a perfect gas with constant specific heats and Prandtl number through a 
two-dimensional cascade, as shown in Fig. 1. In the following discussion, all flow variables 
and spatial coordinates are dimensionless. Lengths have been scaled with respect to the 
blade chord [L*), density, velocity and viscosity with respect to their inlet freestream values 
[p*_ , V*^ and respectively), pressure with respect to twice the inlet freestream dynamic 
pressure (p* ^ V*^), and temperature with respect to the square of the inlet freestream speed 
divided by the specific heat at constant pressure (V_*^/c*). Here the superscript * denotes 
a dimensional quantity and the subscript — oo refers to the prescribed freestream conditions 
far upstream. 
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3. Inviscid Region 


Since the inviscid flow is assumed to be isentropic and irrotational, a velocity potential, 
$, exists and is governed by the field equation 

A 2 V 2 <t> = V(V$) 2 /2 . (3.1) 

The speed of sound propagation, A, the fluid pressure, P, density, p, and temperature, T, 
are related through Bernoulli’s equation and the isentropic relations as follows: 

= ( 7 M’_ Pp-'V-’ = p -' > = = l-b-=4>M?J(V*) 2 -l] , (3.2) 

where M is the Mach number and 7 is the specific heat ratio of the fluid. The inviscid flow 
is determined as a solution of Eq. (3.1) subject to a flow tangency condition at each blade 
surface, cascade periodicity conditions upstream and downstream of the blade row, jump 
conditions on normal velocity and pressure across blade wakes, and the appropriate uniform 
flow conditions far upstream of the blade row. Usually, a Kutta condition is applied at blade 
trailing edges in lieu of specifying the exit flow angle. Finally, far downstream of the blade 
row global mass conservation is enforced, accounting for blockage effects due to the viscous 
layers. For the flows considered here, the inlet and exit velocities are subsonic. 

The specific forms of the blade and wake conditions follow from an asymptotic matching 
of the outer inviscid and the inner viscous-layer equations [1]. Thus, the inviscid solution for 
the normal velocity at a blade surface must match the viscous solution for this velocity at 
the outer edge of the viscous layer. It follows, after carrying out the asymptotic matching, 
that 

V$ • n| s = p~ 1 d(p e u e S)/ds , (3.3) 

where S denotes a reference blade surface (see Fig. 1), p e and u e are the inviscid density and 
velocity at this surface (or, the viscous density and streamwise velocity component at the 
edge of the viscous layer) and <5 is the boundary-layer displacement thickness. The quantities 
s and n denote the arc distance along the blade (positive in the downstream direction and 
zero at the leading-edge stagnation point) and the local unit normal vector directed outward 
from the surface, respectively. 

Two types of terms arise from the wake matching conditions [1], one due to the displace- 
ment thickness effect of the wake and the other due to the wake curvature effect. The first 
leads to the requirement that the inviscid solution for the normal component of velocity 
must be discontinuous with jump given by 

[V$l • n| w = p~ 1 d(p e u e 6 w )/ds , (3.4) 

where [ ] denotes the difference in a quantity (upper minus lower) across the wake, n is the 
upward pointing unit normal vector to the reference wake streamline (i.e., W in Fig. 1), and 
6 W is the displacement thickness of the complete wake. The wake curvature effect gives rise 
to a static pressure difference across the wake. The requirement that the outer inviscid flow 
match this pressure difference leads to the condition 

[P] w = K PeU 2 e {6 w + 6 W ) , (3.5) 
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where 6 W is the momentum thickness of the complete wake and k is the curvature of the 
wake which is taken as positive when the reference wake streamline is concave upwards. 

A complication arises in that the location of the reference wake streamline is unknown 
a priori ; however, to within lowest order, the wake conditions can be referenced to any 
arbitrary curve emanating from the trailing edge and lying within the actual viscous wake 
[8]. In the present study, the reference wake streamline is taken to be the aft stagnation 
streamline as determined from the pure inviscid solution. This is adequate except in extreme 
cases where the location of the stagnation streamline is significantly altered by viscous effects. 
In this case, it is possible to periodically update the location of the wake streamline during 
the calculation, although this has not been done in the present study. 

The boundary- value problem posed by Eqs. (3.1) and (3.2) is solved using Newton iter- 
ation. This is accomplished by substituting for the potential function $ in Eqs. (3.1) and 
(3.2) using the expression 

$ = $n + <^n , (3-6) 

where n represents the Newton iteration count, $ n is the initial guess or previous iteration 
value for the potential function at each point in the field and 4> n is the correction to $ n . The 
variable <f> n is the quantity solved for during each Newton iteration. Higher-order terms (i.e., 
products of < j> n ) are dropped, resulting in a linear equation for 4> n . This equation is discretized 
and the resulting matrix system of algebraic difference equations is inverted directly using 
lower-upper decomposition and Gaussian elimination. The Newton iterations are continued 
until ||^ n || < e, where e is a user specified convergence criterion (e <C 1). The SFLOW 
analysis is described in detail in Ref. [17]. 
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4. Viscous Layer 


The flow in the inner or viscous region is assumed to be governed by Prandtl’s viscous- 
layer equations. After introducing the scaled normal coordinate, n, and normal velocity 
component, v, where 

h = Re l t 2 n and v = Re l ^ 2 v , (4-1) 

the continuity and streamwise momentum equations have the form 


d(pu) + d{pv) 


ds 


dn 


= 0 


and 


. du ,dv, du e d du 

p[u a; + "M 1 - '’*“‘17 = 


(4.2) 


(4.3) 


where u is aligned with the local surface tangent and u > 0 when the flow is in the direction 
of increasing s. Here we assume that the flow in the viscous layer is adiabatic at unit Prandtl 
number; thus, the energy equation reduces to the requirement that the total enthalpy of the 
fluid, H = T + ^ 2 /2, must be constant across the viscous layer. In Eq. (4.3) the subscript e 
refers to fluid properties at the edge of the viscous layer, and the effective turbulent viscosity, 
[i T , is defined by 

V-t = /* + 7r e . ( 4 - 4 ) 

where e is the turbulent eddy viscosity, 7 r is the longitudinal intermittency factor and p 
is the molecular viscosity, which is assumed to be a function of temperature alone. The 
eddy-viscosity model employed in the present study for blade surface boundary layers is 
that of Cebeci and Smith [20], modified to account for separated flow [7]; in the wake, the 
model of Chang, et al. [21] is used. In the present study, the location at which instantaneous 
transition occurs is specified. 

The foregoing field equations govern the flow in the viscous layers along the upper and 
lower surfaces of the blades and in the blade wakes. They are solved subject to conditions at 
the edges of the viscous layer, on the airfoil surface, and along the reference wake streamline, 
i.e., 

u — ► u e for h — ► oo , s > 0 , (4-5) 

u = v — 0 for h = 0, 0 < s < (4-6) 

and 


v = 0 for h = 0, s > s 


TE 


(4.7) 

respectively, where s± E are the trailing-edge values of the upper- (+) and lower-surface ( — ) 
arc-length coordinates measured from the leading-edge stagnation point. The condition 
expressed by Eq. (4.5) is also applied along the wake streamline for h — ► — oo. Equations 
(4.6) and (4.7) imply that the curve n = 0 corresponds to the blade surfaces and reference 
wake streamlines, respectively. 

The displacement and momentum thicknesses of the viscous layers (S and 9 , respectively) 
are needed to determine the effect of viscous displacement and wake curvature on the outer 
inviscid flow (Eqs. (3.3)-(3.5)). They are defined by 

pu 


S(s) = Re~ 1/2 [ C 
Jo 


(1 - -!—)dh 


Pe^e 


(4.8) 
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and 


(4.9) 


Jo Pe^e Ue 


where the zero lower bound on the integrals is replaced with —oo in the wake. 

The independent and dependent variables appearing in the viscous-layer equations are 
transformed using a modified form [7] of the Levy-Lees transformation [22]. Thus, the 
independent variables are given by 



Pe^ePe^ds 


and 


Tj — Pe^e 


(20 


- 1/2 


J p ! pt 

Jo 


dh 


and the dependent variables are defined by 


F = u/u e and / = ( 2 £) . 


(4.10) 


(4.11) 


The quantity e is the value of the “outer” eddy-viscosity parameter (Eq. (4.4)), at each 
s-station, as defined in the Cebeci-Smith model [20]. It appears in the definition of £ in 
order to maintain a nearly constant 77 value for the edge of the boundary layer in a turbulent 
flow [23]. The Levy-Lees transformation permits the leading-edge stagnation point similarity 
solution to be recovered and reduces the truncation error of the viscous solution over that 
associated with the use of primitive variables. The quantity if; is the compressible stream 
function defined in terms of the surface and wake streamline tangential and scaled normal 
velocity components u and u, respectively, from the relations pu = dip/dh and pv = —dxl>/ds. 

The continuity and momentum equations transform to: 

(4.12) 



and 



(4.13) 


Here t. = pp./ p e p e , £ is the ratio of the turbulent to molecular viscosity coefficients and /? is 
the pressure gradient parameter, defined by 


2 £ du e 

u e d£ 


(4.14) 


The quantity 0 is equal to T/T e , the local to edge static temperature ratio, which is related 
to F as follows: 

e = 1 + ^=J-M e 2 (l - F 2 ) . (4.15) 

The molecular viscosity coefficient, //, is determined using the Sutherland viscosity law; i.e., 


P 

P — oo 


T \ 3/2 T_^ + Tc 
T. oo) T + T c 


(4.16) 


Here p_oo is the molecular viscosity at the temperature T - and Tc is a constant, which for 
air has a dimensional value of 110°K [24]. 
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The following boundary conditions are applied: the blade and wake reference surfaces 
(at 77 = 0) are both streamlines, therefore the stream function is a constant, so that, without 
loss of generality, we can set 

/ = 0 at tj = 0 . (4.17) 

On the blade surface the no-slip condition applies so that 

F = 0 at 77 = 0, f < t TE , (4.18) 

where &E is the trailing-edge value of £. At the edges of the surface and wake viscous layers 
the remaining boundary conditions are applied, namely 


F(vt) = F^;) = 1 , (4.19) 

where 77+ and tj~ are the upper and lower viscous-layer edge values of the wake 77-coordinate, 
respectively; the ± superscripts are not needed for the surface boundary layers since each 
has only one edge. The boundary condition given by Eq. (4.19) enforces the approach of the 
flow variables to their correct edge values. Equation (4.19) is consistent with the assumption 
that the static pressure jump across the wake is negligible. 

The viscous-layer equations are parabolic in the ^-direction and therefore require initial 
conditions. These are provided by a similarity solution obtained at the leading-edge stag- 
nation point, £ = 0, by solving the above equations with 0=1. The blade surface and 
wake solutions are obtained by space-marching in the downstream direction. As discussed 
in §2, a complete IVI calculation requires the ability to solve the viscous-layer equations in 
both the “direct” mode, where the pressure gradient parameter 0 is specified and the dis- 
placement thickness 8 is unknown, and in the “inverse” mode, where 8 is specified and 0 is 
unknown. In particular, an inverse calculation is required to predict viscous-layer separation. 
The equations are solved in the direct mode near the leading edge, and in the inverse mode 
downstream of an axial station chosen to ensure that the inverse mode is initiated upstream 
of a separation point. The wake is calculated using the inverse mode. 

In the direct viscous-layer calculation, the value of 0 is determined by the inviscid anal- 
ysis and the displacement thickness is obtained from the viscous analysis. In the inverse 
procedure, the value of the mass deficit parameter, m = p e u e 8, is specified and the edge 
values of the variables u e , M e , etc. are obtained as part of the viscous-layer solution. An 
expression for / at the edge of the viscous layer is derived by integrating Eq. (4.12) across 
the boundary layer and employing the definition of 8 (Eq. (4.8)), i.e., 


/( 7?e) = h + 77e - 


m 


>/sr 


(4.20) 


where 


h = r\e-i)dri . 

Jo 


(4.21) 


Equation (4.20) is used on blade surfaces to impose the specified value of m through the 
corresponding value of / at the outer edge of the viscous mesh, 77 = T] e . The parameter 
h is specified, and is lagged from the previous global iteration. Similarly, in the wake an 
expression for the jump in the stream function between the viscous-layer upper and lower 
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edges is obtained by integrating Eq. (4.12) across the entire wake (from r/ e to 77 +) and using 
the definition of the displacement thickness, i.e., 

777 

f(vt ) - f(v7) = + vt - V7 - > (4-22) 

where 

rvi a 

/i W = / [6 — l)dr) 

Equations (4.22) and (4.23) are used in the same way that the corresponding equations were 
used on the surfaces to impose m w = (/? e u e <5) w and hw. In the inverse mode the quantities ft 
and u e are unknown. Thus, a supplemental equation relating these two variables is needed. 
This relation is obtained by discretizing Eq. (4.14), which defines ft in terms of u e . 

The discretized governing equations, boundary and auxiliary conditions, Eqs. (4.12)- 
(4.23), are quasi-linearized and the resulting coupled tridiagonal system of algebraic equa- 
tions is solved at each s-station, using a fixed-point iteration to update the nonlinear terms. 
The inversion algorithm used in the wake is modified to account for the application of one 
boundary condition (Eq. (4.17)) at rj = 0 and the others (Eq. (4.19)) at the upper and lower 
edges of the viscous layer, as well as to account for the application of a jump condition on 
/ (Eq. (4.22)) between the upper and lower edges. Finally, the so-called FLARE approxi- 
mation, which prevents instabilities in the viscous-layer solution due to axial flow reversal, 
is applied by turning off all of the convective terms in the momentum equation wherever 
F < 0. The details of the viscous-layer numerical analysis are provided in Appendix A of 
this report. 


(4.23) 
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5. Inviscid/ Viscid Interaction 


The IVI approach used here determines the complete flow field by iteratively updating 
the mass deficit parameter, m(s), which affects the inviscid and viscous solutions through 
their respective boundary conditions. For an arbitrary m distribution, two different surface 
and wake streamline velocity distributions generally result — one, u ej (s ), from the inviscid 
calculation and one, u ev (s), from the viscous-layer calculation. The objective is to deter- 
mine a converged inviscid/viscid interaction solution by finding the mass deficit parameter 
distribution that minimizes the difference between the u ei and u ev distributions. In this in- 
vestigation the so-called “semi-inverse” iteration procedure of Carter [25] is used to update m 
at every streamwise mesh station on the blade and wake surfaces. The term “semi-inverse” 
refers to the use of an inverse method (i.e., displacement thickness specified and pressure 
determined) to solve the viscous-layer equations coupled to a direct method (i.e., geometry 
specified and pressure determined) to solve the inviscid equations. We set 

m n+1 = rn*[\ + - 1)] , (5.1) 

where the superscript n is the global iteration count and a; is a relaxation parameter. The 
solution is considered to be converged when 

max |u ev , - u e , \/u e . < e , i = 1, ...,IE , (5.2) 

1 * 11 

where the value of t is specified by the user and IE is the number of streamwise mesh 
stations. Equation (5.2) is applied on both blade surfaces and along the wake. The viscous- 
layer solution is obtained at the locations corresponding to the intersections of the inviscid 
mesh with the blade and wake streamline surfaces, which avoids the need for interpolation. 

Because a major objective of this study is to develop an efficient analysis, various tech- 
niques for accelerating convergence were examined. We found that one of the most effective 
approaches for reducing the CPU time needed to obtain a converged IVI solution is to use 
the largest value of the relaxation parameter ui for which the iterative procedure remains 
stable. It was observed during this study that the inviscid surface velocity distribution, u en 
undergoes relatively little change from the initial purely inviscid solution to the final con- 
verged IVI solution. This is in contrast to the significantly larger change between the initial 
and final viscous-layer edge distribution, u ev . This observation prompted the introduction of 
a subiteration loop in which the viscous equations are solved repeatedly (“A^^” times) during 
each global iteration. Thus, Eq. (5.1) is applied N v times during a single global iteration 
with u ej frozen at its most recent value while u ev is recalculated during each subiteration 
by solving the viscous-layer equations using the latest m distribution. The value of N v is 
a user-specified input; the standard iteration procedure is recovered when N v = 1. This 
strategy is only effective in reducing the total CPU time if the number of global iterations 
needed to obtain a converged IVI solution, N g , can be reduced enough to more than balance 
the increased computational effort needed for the additional viscous calculations performed 
during each global iteration. Of the three cascade configurations examined in §6, only the 
turbine cascade benefitted from this approach, as discussed in the next section. 

The semi-inverse iteration procedure is illustrated in Fig. 2, where the dashed line corre- 
sponds to the viscous subiteration technique described above. 
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6. Numerical Examples 


The foregoing inviscid/viscid interaction analysis has been applied to the following cas- 
cade configurations: a compressor exit guide vane (EGV), a high-speed compressor (HSC) 
and a turbine cascade. For the compressor cascades the behavior of the surface and wake 
pressure coefficient, Cp = (P — F_ 00 )/2 will be examined. The surface Mach number dis- 
tribution will be examined for the turbine cascade. The displacement thickness, 8, and 
the surface shear stress, t w = Rt~ x [idujdn | n=0 , distributions will be presented for all three 
configurations. The IVI solutions for the compressor cascades will be evaluated through 
comparisons with Navier-Stokes solutions; that for the turbine with experimental measure- 
ments. In addition, the predicted values of the total pressure loss and exit flow angle are 
presented for the EGV cascade. These were obtained using the mixing analysis of Stewart 
[28]. Alternatively, total pressure loss and exit flow angle could be based on the predicted 
average flow conditions at the exit plane of the computational domain. Both methods have 
been implemented in SFLOW-IVI and give results that are in very good agreement for the 
flows considered in this study. Finally, the performance of the SFLOW-IVI analysis, i.e., its 
efficiency and convergence properties, will also be discussed. 

In all of the calculations described here, the SFLOW-IVI analysis was applied using a 
convergence tolerance, e, of 0.001 (see Eq. (5.2)). The inviscid meshes that were used for both 
compressor cascades have the same dimensions, i.e., there are 90 axial and 31 circumferential 
lines, with 24 axial lines upstream of the leading edge, 41 lines intersecting the blade surfaces 
and 25 lines aft of the trailing edge. The viscous-layer analysis employed a total of 81 blade 
and 25 wake streamwise grid lines, with 71 grid lines across the surface boundary layer and 
141 grid lines across the wake. The inviscid mesh used for the turbine cascade had 150 axial 
and 31 circumferential lines, with 39 points upstream of the leading edge, 51 points along 
each blade surface and 60 points along the wake. A total of 101 surface and 25 wake stations 
were used in the viscous-layer analysis. The normal mesh had the same dimensions as those 
used for the compressor cascades. For the cases considered in this study, the wake curvature 
effect was regarded to be negligible; thus, [PJ^, was set equal to zero, cf Eq. (3.5). 

The inviscid solutions were obtained on a “streamline” type H-mesh, rather than on the 
“sheared” H-mesh described in Refs. [14] and [17]. The SFLOW analysis of Ref. [17] was 
modified to use the streamline H-mesh developed by Hall and Verdon [26]. Thus, prior to 
initiating an IVI calculation, an inviscid solution is obtained on a sheared H-mesh. This 
solution is then used to generate a streamline H-mesh, in which one family of mesh lines 
corresponds to the streamlines of the inviscid flow solution obtained using the simple sheared 
H-type mesh, and the second family consists of lines which are “nearly” orthogonal to the 
first set. The principal advantage of the streamline H-mesh over the sheared mesh is an 
improved resolution of the flow in the blade leading-edge regions. 

An alternative to this procedure is available in SFLOW and was used for one of the 
cases described below, i.e., the turbine cascade. In the turbine case, a useful streamline 
H-mesh could not be produced from the solution obtained using a sheared H-mesh, because 
the latter provides inadequate resolution of the leading-edge region. To remedy this, the 
initial inviscid solution was obtained on a “composite” mesh constructed by overlaying a 
surface fitted C-mesh, generated locally in the leading-edge region, over the global sheared 
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H-mesh. A detailed description of this procedure, as applied to linearized unsteady flows, 
can be found in Ref. [27]. 

6.1 Compressor Exit Guide Vane 

The EGV cascade consists of 12 percent thick, highly cambered, modified NACA airfoils 
[26]. It has a stagger angle, 0, of 15 deg, a gap-chord ratio, G, of 0.6 and operates at a pre- 
scribed inlet Mach number, M_oo, and inlet flow angle, H-oo, of 0.3 and 40 deg, respectively. 
Calculations were performed for an inviscid flow, and for viscous flows at Reynolds numbers 
of 10 5 and 10 6 . Instantaneous transition from laminar to turbulent flow was assumed to 
occur at one percent of the arc distance measured from the leading-edge stagnation point to 
the trailing edge on both the suction and pressure surfaces of the blades. This is a physically 
realistic assumption for most high Reynolds number flows and is consistent with the usual 
practice in high Re Navier-Stokes calculations in which instantaneous transition is assumed 
to occur at blade leading edges. A streamline mesh is depicted in Fig. 3, where three adja- 
cent EGV blade passages are shown. For the purpose of illustration, the mesh shown in this 
figure has approximately half as many axial and circumferential grid lines than were used 
for the actual calculations. 

Results of the inviscid and IVI calculations are shown in Fig. 4. The blade and wake 
pressure and displacement thickness distributions are shown in Figs. 4(a) and 4(b), respec- 
tively, and the surface shear-stress distributions along the blade in Fig. 4(c). The expected 
approach of the viscous to the inviscid solution as Re is increased is evident in Fig. 4(a). 
The rate of growth of the suction-surface displacement thickness increases with increasing x 
as the viscous-layer separation point is approached. In the wake the half-wake displacement 
thickness, <*>vv/ 2 is plotted in Fig. 4(b). As expected, its value just aft of the trailing edge is 
approximately equal to the mean of the upper- and lower-surface trailing-edge displacement 
thicknesses. As shown in Fig. 4(c), a suction-surface separation bubble (r w < 0) exists and 
spans approximately 14 percent of chord at Re = 10 6 , and about 24 percent of chord for 
Re = 10 5 . The decrease in the extent of the separation bubble as Re is increased is consistent 
with the behavior expected for turbulent flows. 

The surface pressure, displacement thickness and surface shear-stress distributions pre- 
dicted by SFLOW-IVI are compared in Fig. 5 with results obtained using the Navier-Stokes 
analysis of Dorney, et al. [29] for the Re = 10 6 case. The N-S analysis uses the Baldwin- 
Lomax turbulence model [30], which is very similar to the Cebeci-Smith model used in 
SFLOW-IVI. Good agreement between the results of the two procedures is obtained over 
most of the blade surface for all three quantities. However, the agreement deteriorates in the 
vicinity of the trailing edge. This is caused by the use of an O-mesh around the blades in 
the Navier-Stokes analysis. The O-mesh topology is not well-suited to wedge-shaped trailing 
edges like those found in the EGV cascade. Both analyses predict separation (t w < 0) near 
the trailing edge, and give almost identical predictions for the location of the separation 
point (t w = 0; see Fig. 5(c)). 

To test the robustness of the SFLOW-IVI analysis, additional calculations were carried 
out for M_oo = 0.3, Re = 10 6 and a wide range of inlet flow angles, 36 deg < <*> < 54 deg. 

The transition point locations were held fixed at s/sje = 0.01 for all values of fl-oo. This 
location is the same as that reported above for the baseline (ft_oo = 40 deg) calculation. 
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The results are shown in Fig. 6, where the predicted total pressure loss parameter, u7 = 
(Pt_ x — Pt+^/iPt-oo — P- oo) (where P t is the total pressure), exit flow angle, fl +00 , and 
separation point location, x sep , are plotted as functions of fi_oo- Solutions could not be 
obtained for fl_ <*, < 36 deg because a small supersonic region, which could not be treated 
using the present version of SFLOW-IVI, formed near the leading edge of each blade. At 
Cl-oo = 54 deg, the viscous layer is approaching stall, with the separation region spanning 
approximately 35 percent of chord. Above 54 deg the solution would not converge due to 
a numerical instability. This is consistent with the known stability properties of the semi- 
inverse IVI iteration procedure when applied to flows with large-scale separations [31]. 

The total pressure loss parameter and the exit flow angle are plotted versus ^ in 
Figs. 6(a) and 6(b), respectively. There is a wide range of inlet flow angles over which the 
loss remains relatively low, while ZJ increases rapidly as the inlet flow angle is increased above 
50 deg. The latter corresponds to the rapid inflation of the separation region with increasing 
fi-oo for f!_oo > 50 deg, which can be seen in Fig. 6(c). A striking similarity exists between 
the variations in — fl +oc , and x sep as Cl_ 00 is varied, as is apparent from comparing the results 
shown in Figs. 6(b) and 6(c). The streamwise growth of the separation bubble as is 

increased is accompanied by a similar increase in the suction-surface displacement thickness 
in the vicinity of the trailing edge. This produces a thickened displacement body (i.e., blade 
plus displacement thickness), reducing the effective camber of the blade and thus, the loading 
it produces. As a direct consequence there is a reduction in the turning of the flow, i.e., an 
increase in fl +00 . 

The predicted streamline patterns indicating the size of the trailing-edge separation bub- 
ble for fi_oo = 36 deg, 45 deg and 54 deg are shown in Fig. 7. The separation bubble 
grows slowly in the range 36 < fl-oo < 45 deg, and much more rapidly between 45 and 54 
deg. The “decambering” effect produced by the growth of the separation bubble is clearly 
illustrated by these results. The kinks that appear in the streamlines near the trailing edge, 
depicted in Fig. 7, require some explanation. Since the blade trailing edge is wedge-shaped, 
the surface coordinate line formed by the blade surface and reference wake streamline has a 
geometric singularity or “kink” at the trailing edge. This singularity influences the solution 
throughout the trailing-edge region as shown in the streamline plots in Fig. 7. Because this 
singular behavior is highly localized, its effect on the overall flow field solution is negligible. 
An additional factor contributing to the kinks is the use of a relatively coarse streamwise 
mesh aft of the trailing edge. This results in an abrupt change in the flow variables across 
the trailing-edge point. 

6.2 High-Speed Compressor 

The HSC cascade consists of cambered, modified NACA 0006 airfoils [32]. This cascade 
operates at high-subsonic inlet conditions, i.e., M_ <*, = 0.7 and O_oo = 55 deg, and has 
a blade spacing and a stagger angle of unity and 45 deg, respectively. As was done for 
the EGV cascade, two values of the Reynolds number have been considered, 10 5 and 10 6 . 
Instantaneous transition was assumed to occur at ten and at one percent of the surface arc 
length for Re = 10 5 and 1 0 6 , respectively, on both the suction and pressure surfaces of each 
blade. The HSC cascade solution for Re = 10 5 was found to be sensitive to the specified 
location of transition. The solution could not be converged with transition specified at 
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one percent of arc length, however, specifying transition farther downstream permitted the 
solution to converge. The reason for this sensitivity has not yet been determined and needs 
to be examined in future work. A mesh which has, for the sake of clarity, a lower grid point 
density than was used for the actual calculations is shown in Fig. 8. Inviscid and viscous 
calculations were performed and the results are presented in Figs. 9 and 10. 

The predicted pressure and displacement thickness distributions along the blade surfaces 
and the wake are shown in Figs. 9(a) and 9(b), respectively. The behavior of both quantities 
is similar to that observed for the EGV. The surface shear-stress distributions shown in 
Fig. 9(c) indicate that the extents of the suction-surface separation bubbles are smaller than 
those predicted for the EGV cascade, decreasing from approximately 20 percent to about 
8 percent of chord as the Reynolds number is increased from 10 5 to 10 6 . The kinks in the 
shear-stress distributions for the Re = 10 5 case are associated with transition. 

The surface pressure coefficient, displacement thickness and shear-stress distributions 
obtained for Re = 10 6 using SFLOW-IVI are compared in Fig. 10 with those obtained using 
the Navier-Stokes analysis of Ref. [29]. The agreement is excellent except in the immediate 
vicinity of the trailing edge. The two analyses give almost identical predictions for the 
location of the separation point. Again, the differences in the two solutions are attributed 
to the use of an O-mesh for wedge-shaped trailing-edge geometries in the Navier-Stokes 
analysis. 

6.3 Turbine 

The turbine cascade considered here is the Fourth Standard Configuration described in 
the study of Fransson and Suter [33]. The blade geometry is shown in Fig. 11 and was 
obtained by modifying the original blunt trailing-edge geometry to produce a wedge-shaped 
trailing edge while retaining the original chord length, as discussed in Ref. [17]. As for the 
cases discussed above, the mesh shown in the figure has fewer grid lines than were used in the 
actual calculation. The streamline mesh employed for the turbine calculation was obtained 
from an inviscid solution calculated using the composite mesh analysis discussed in §4 and 
in Ref [27], 

The blade spacing and stagger angle for the turbine cascade are 0.76 and 56.6 deg, 
respectively, and the inlet Mach number and flow angle are 0.205 and 45 deg, respectively. 
The value of A/_ <*> has been adjusted from the experimentally measured value of 0.190 to 
improve the agreement with the measured pressure distribution. The calculation was carried 
out at a Reynolds number of 5 x 10 5 with instantaneous transition occurring at 10 percent of 
the surface arc length along both surfaces. A converged solution could not be obtained for 
the turbine cascade if the location of transition was specified to be too close to the leading 
edge. This was consistent with the behavior observed for the HSC cascade calculation at 
Re = 10 5 . 

The IVI solution was obtained in 12 global inviscid/viscid iterations. The viscous subit- 
eration procedure described in §5 was very effective for this case, reducing the CPU time 
needed to converge the calculation from 1371 seconds without subiteration (requiring 115 
global iterations) to 224 seconds (in 12 global iterations), using four viscous subiterations 
(i.e., N v = 4) during each global IVI iteration. 

The computed and measured blade surface Mach number distributions are shown in 
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Table 1. Summary of SFLOW-IVI CPU times, tc, for different cascade configurations. 


Configuration 

LJ 

N v 

N g 

tc (secs) 

EGV, Re = 10 6 

1.20 

1 

24 

197 

EGV, Re = 10 s 

0.85 

1 

38 

277 

HSC, Re = 10 6 

1.20 

1 

27 

203 

HSC, Re = 10 5 

0.80 

1 

40 

296 

Turbine 

0.55 

4 

12 

224 


Fig. 12(a). Viscous effects produce a nearly uniform decrease in the suction-surface Mach 
number distribution aft of x fa 0.4, while the pressure-surface distribution is almost unaf- 
fected. The agreement between the IVI solution and the experimental data is reasonable — 
the disagreement in the trailing-edge region may be partly attributable to the local geom- 
etry modification mentioned above. It is difficult to draw definitive conclusions regarding 
the experimental comparison because the solution for this case is particularly sensitive to 
the inviscid mesh. The predicted displacement thickess and surface skin-friction coefficient 
distributions are shown in Figs. 12(b) and 12(c), respectively. No separation was predicted 
for the conditions considered here, although the suction-surface viscous-layer is close to sep- 
aration at the trailing edge. 

6.4 Timing Study and Convergence Behavior 

Because the development of an efficient analysis has been a major objective of this an- 
alytical effort, a timing study was conducted for the three cascade configurations examined 
here. This provides both a measure of the computational effort presently required to obtain 
solutions using SFLOW-IVI, and benchmarks against which future efforts to improve effi- 
ciency can be compared. The results are summarized in Table 1. In addition to the CPU 
time ( tc), the relaxation factor (u;), the number of viscous subiterations used ( N v ) and the 
number of global iterations (N g ) required to converge the IVI solution using a tolerance 
level, e, of 0.001, are given in Table 1. The execution times were determined using the nearly 
optimal value of u>, which was obtained by trial and error. 

The calculations were carried out on an HP-Apollo 720 workstation where SFLOW-IVI 
has been compiled using an optimizing preprocessor. No attempt has been made to “tune” 
the code to take advantage of special features of the optimizer. The times given in Table 
1 are CPU times for the portion of the calculation associated with the IVI iteration loop. 
Any overhead associated with initialization of the data structure, generation of the mesh and 
calculation of the initial inviscid solution is not included. However, this overhead amounts 
to a small percentage of the overall CPU time required by the present analysis. Note that 
each of the solutions was obtained in less than five minutes. 

The convergence behavior of two parameters of interest to compressor blade designers 
has been examined here to determine if the solutions discussed above are sufficiently con- 


17 


verged, and whether a different measure of convergence than that given by Eq. (5.2) would 
be more appropriate. For the two compressor cascades, the total pressure loss parameter 
U and exit flow angle fI +00 were monitored during the IVI iterations. We have found that, 
for most engineering purposes, the solutions obtained herein could be considered converged 
at a significantly lower iteration count than was needed to satisfy the convergence criterion 
(e = 0.001). Thus, even greater efficiency could be achieved in many cases by measur- 
ing convergence by the degree to which the parameters of interest have approached their 
“asymptotic” values. This is demonstrated by the results presented in Fig. 13, which show 
the behavior of ZJ and 17 +00 , respectively, as functions of the iteration count for the EGV 
cascade operating at Re = 10 6 and ff-oo = 40 deg. This behavior is typical of that observed 
for all of the cases studied herein. The solution for the case illustrated in Fig. 13 converged to 
e = 0.001 within 24 iterations while the asymptotic value, indicated by the dashed horizontal 
line, was determined by converging the solution to t = 0.0001, for which 41 iterations were 
needed. For engineering purposes, this solution could be considered to be converged after 
about 15 iterations, for which tc is approximately 120 seconds. 

The CPU times required to obtain the present IVI solutions are significantly lower than 
those needed to obtain Navier-Stokes solutions of comparable accuracy. The latter currently 
require approximately one to three hours of CPU time on a modern workstation. However, 
to achieve such CPU times considerably fewer grid points (typically between 8 and 15) are 
used across the surface boundary layers than are used in IVI calculations. For example, 71 
points were used for the SFLOW-IVI calculations presented in this report. Thus SFLOW-IVI 
requires one to two orders-of-magnitude less CPU time than current Navier-Stokes analyses 
to achieve viscous solutions for 2-D cascade flows. 
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7. Concluding Remarks 


In the present study, existing nonlinear inviscid and inverse viscous-layer analyses have 
been extended and coupled to provide a strong inviscid/viscid interaction solution capability 
for two-dimensional cascade flows. This IVI solution procedure can be used to predict the 
effects of local strong interactions, including trailing-edge/near-wake interactions and small- 
to moderate-scale viscous-layer separations, on cascade performance. The present analysis 
is restricted to subsonic flows, but it can be extended to treat transonic flows. 

The SFLOW-IVI analysis has proven to be accurate, efficient and robust. During this 
study the analysis was applied to compressor exit guide vane, high-speed compressor and 
turbine cascades. The analysis has been shown to be able to predict both the detailed 
features of the flow field as well as global quantities such as loss and turning. Very good 
agreement with Navier-Stokes solutions was demonstrated for two of the cases considered 
herein. 

Converged solutions for each of the foregoing configurations examined here were obtained 
in less than five minutes on an HP-Apollo 720 Workstation. It was shown that even lower 
CPU times could be obtained by basing convergence on the global quantities of interest to 
the engine designer. The robustness of the SFLOW-IVI analysis was demonstrated through 
application to a wide range of inlet conditions, including cases in which large-scale separa- 
tions, spanning up to 35 percent of chord, occurred. It should be noted that CPU times for 
the most severe cases were on the order of 15 to 20 minutes. 

A number of issues still need to be addressed in order to improve the accuracy of the 
present steady IVI analysis, and to expand its range of applicability. Among them are the 
inclusion of quasi-three-dimensional effects (i.e., streamtube contraction and radius change), 
the incorporation of predictive models for determining the transition from laminar to turbu- 
lent flow and the addition of a procedure for updating the location of the wake streamline 
during the global iteration process. With respect to transition, at present we are specifying 
the transition location, but for a truly predictive viscous calculation (either IVI or Navier- 
Stokes) the transition location should be determined as part of the complete solution. In 
addition, the overall utility of the SFLOW-IVI analysis for design-system applications needs 
to be explored through further testing and validation. Finally, as this effort continues, the 
focus will increasingly turn towards the development of an unsteady strong inviscid/viscid 
interaction capability based, as much as possible, on the complementary SFLOW-IVI and 
LINFLO analyses. 


References 


1. Melnik, R. E., “Turbulent Interactions on Airfoils at Transonic Speeds — Recent Devel- 
opments,” AGARD-CP-291, Chapter 10, 1980. 

2. Lock, R. C. and Firmin, M. C. P., “Survey of Techniques for Estimating Viscous Effects 
in External Aerodynamics,” Royal Aircraft Establishment Technical Memorandum Aero 
1900, April 1981. 

3. Girodroux-Lavigne, P. and Le Balleur, J. C., “Unsteady Viscous-Inviscid Interaction 
Method and Computation of Buffeting over Airfoils,” ONERA T.P. No. 1987-58, 1987. 

4. Howlett, J. T., “Efficient Self-Consistent Viscous-Inviscid Solutions for Unsteady Tran- 
sonic Flow, AIAA Paper 85-0482, 1985. 

5. Caspar, J. R., “Unconditionally Stable Calculation of Transonic Potential Flow through 
Cascades Using an Adaptive Mesh for Shock Capture,” ASME Journal of Engineering 
for Power , Vol. 105, No. 3, 1983, pp. 504-513. 

6. Vatsa, V. N. and Verdon, J. M., “Viscous/Inviscid Interaction Analysis of Separated 
Trailing-Edge Flows,” AIAA Journal , Vol. 23, No. 4, 1985, pp. 481-489. 

7. Barnett, M. and Verdon, J. M., “Viscid/Inviscid Interaction Analysis of Subsonic Turbu- 
lent Trailing-Edge Flows,” AIAA Journal , Vol. 25, No. 9, September 1987, pp. 1184-1193. 

8. Veldman, A. E. P., “The Calculation of Incompressible Boundary Layers with Strong 
Viscous-Inviscid Interaction,” AGARD-CP-291, Chap. 12, 1980. 

9. Hansen, E. C., Serovy, G. K. and Sockol, P. M., “Axial-Flow Compressor Turning Angle 
and Loss by Inviscid-Viscous Interaction Blade-to- Blade Computation,” ASME Paper 
No. 79-GT-5, 1979. 

10. Janssens, P. and Hirsch, Ch., “A Viscid Inviscid Interaction Procedure for Two-Dimen- 
sional Cascades,” AGARD-CP-351, Chapter 3, 1983. 

11. Calvert, W. J. and Herbert, M. V., “An Inviscid-Viscous Interaction Method to Predict 
the Blade-to-Blade Performance of Axial Compressors,” Aeronautical Quarterly, Vol. 
XXXI, Part 3, 1980, pp. 173-196. 

12. Barnett, M. and Verdon, J. M., “Theoretical Prediction of High Reynolds Number Vis- 
cid/Inviscid Interaction Phenomena in Cascades,” Proceedings of the Fourth Symposium 
on Numerical and Physical Aspects of Aerodynamic Flows, Long Beach, California, 1989. 

13. Barnett, M., Hobbs, D. E. and Edwards, D. E., “Inviscid-Viscous Interaction Analysis of 
Compressor Cascade Performance,” ASME Journal of Turbomachinery, Vol. 113, No. 4, 
1991, pp. 538-553. 

14. Verdon, J. M. and Caspar, J. R., “A Linearized Unsteady Aerodynamic Analysis for 
Transonic Cascades,” Journal of Fluid Mechanics, Vol. 149, 1984, pp. 403-429. 


20 


15. Power, G. D., Verdon, J. M. and Kousen, K. A., “Analysis of Unsteady Compressible 
Viscous Layers,” ASME Journal of Turbomachinery , Vol. 113, 1991, pp. 644-653. 

16. Barnett, M. and Verdon, J. M., “Analysis of Blade Unsteady Boundary Layers and 
Wakes,” presented at the Sixth International Symposium on Unsteady Aerodynamics, 
Aeroacoustics and Aeroelasticity of Turbomachines and Propellers, Notre Dame, Indiana, 
September 1991; Proceedings to be published by Springer- Verlag. 

17. Hoyniak, D. and Verdon, J. M., “Development of a Steady Potential Solver for Use With 
Linearized, Unsteady Aerodynamic Analyses,” NASA TM 105288, September 1991. 

18. Goldstein, S. “On Laminar Boundary Layer Flow Near a Position of Separation,” Quar- 
terly Journal of Mechanics and Applied Mathematics, Vol. 1, 1948, pp. 43-69. 

19. Stewartson, K., “Is the Singularity at Separation Removable?” Journal of Fluid Mechan- 
ics, Vol. 44, Part 2, 1970, pp. 347-364. 

20. Cebeci, T. and Smith, A. M. 0., Analysis of Turbulent Boundary Layers, Academic 
Press, New York, 1974. 

21. Chang, K. C., Bui, M. N., Cebeci, T. and Whitelaw, J. H., “The Calculation of Turbulent 
Wakes,” AIAA Journal, Vol. 24, No. 2, 1986, pp. 200-201. 

22. Blottner, F. G., “Finite Difference Method of Solution of the Boundary- Layer Equations,” 
AIAA Journal, Vol. 18, No. 2, 1970, pp. 193-205. 

23. Vatsa, V.N., Werle, M.J. and Verdon, J.M., “Viscid/Inviscid Interaction at Laminar and 
Turbulent Symmetric Trailing Edges,” AIAA Paper 82-0165, 1982. 

24. Schlichting, H., Boundary-Layer Theory, Seventh Edition, McGraw-Hill Book Company, 
Inc., New York, 1979, p. 328. 

25. Carter, J. E., “A New Boundary Layer Inviscid Iteration Technique for Separated Flow,” 
AIAA Paper 78-1450, 1978. 

26. Hall, K. C. and Verdon, J. M., “Gust Response Analysis for Cascades Operating in 
Nonuniform Mean Flows,” AIAA Journal, Vol. 29, No. 9, 1991, pp. 1463-1471. 

27. Usab, W. J., Jr. and Verdon, J. M., “Advances in the Numerical Analysis of Linearized 
Unsteady Cascade Flows,” ASME Journal of Turbomachinery, Vol. 113, 1991, pp. 633- 
643. 

28. Stewart, W.L., “Analysis of Two-Dimensional Compressible Flow Loss Characteristics 
Downstream of Turbomachine Blade Rows in Terms of Basic Boundary- Layer Charac- 
teristics,” NACA TN 3515, 1955. 

29. Dorney, D.J., Davis, R.L. and Edwards, D.E., “Investigation of Hot Streak Migration and 
Film Cooling Effects on Heat Transfer in Rotor/Stator Interacting Flows,” UTRC Report 
91-29, April 1992. Prepared under N A VAIR Contract N00140-88-C-0677 (Report 1). 


21 


30. Baldwin, B. S. and Lomax, H., “Thin-Layer Approximation and Algebraic Model for 
Separated Turbulent Flows,” AIAA Paper 78-257, 1978. 

31. Wigton, L.B. and Holt, M., “Viscous-Inviscid Interaction in Transonic Flow,” AIAA 
Paper 81-1003, 1981. 

32. Verdon, J. M. and Hall, K. C., “Development of a Linearized Unsteady Aerodynamic 
Analysis for Cascade Gust Response Predictions,” NASA CR 4308, July 1990. 

33. Fransson, T. H. and Suter, P., “Two-Dimensional and Quasi Three-Dimensional Exper- 
imental Standard Configurations for Aeroelastic Investigations in Turbomachine Cas- 
cades,” Report LTA-TM-83-2, Ecole Polytechnique Federale De Lausanne, Lausanne, 
Switzerland, 1983. 

34. Reyhner, T. A. and Flugge-Lotz, I., “The Interaction of a Shock Wave with a Laminar 
Boundary Layer,” Int. J. Non-Linear Mech., Vol. 8, No. 2, 1968, pp. 173-193. 


22 


List of Symbols 


All physical parameters listed below are dimensionless, as described in §3, General Con- 
cepts. The number(s) in parentheses at the end of each symbol description indicates an 
equation in which the symbol appears. 

Roman 

A Speed of sound propagation, (3.1). 

Cp Pressure coefficient, Figure 4. 

e Turbulent eddy viscosity, (4.4). 

F, f Dependent variables in Levy-Lees transformation, (4.11). 

G Cascade gap-to-chord ratio, Figure 1. 

H Total enthalpy. 

h Viscous-layer integral parameter, (4.21). 

I Viscous-layer parameter, (4.13). 

M Mach number, (3.2). 

m mass deficit parameter, (4.20), (4.22). 

N g Number of global iterations to converge solution, Table 1. 

N v Number of viscous subiterations per global iteration. 

n Scaled viscous-layer normal distance from surface and wake reference streamline, 
(4.1). 

n Unit normal positive when directed outward from blade surface or upward from 
wake reference streamline, (3.3), (3.4). 

P Pressure, (3.2). 

Re Reynolds number. 

S Blade surface, Figure 1. 
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s Surface and wake streamline arc-distance measured from leading-edge stagnation 
point, (4.2, 4. 3). 

T Temperature, (3.2). 

tc CPU time to converge solution, Table 1. 

u, v Viscous-layer velocity components in directions of surface and wake streamline 
tangent and normal, respectively, (4. 1-4.3). 

v Scaled viscous-layer normal velocity component, (4.1). 

W Reference wake streamline, Figure 1. 

x , y Airfoil frame coordinates, Figure 1. 

x sep x-location of separation point, Figure 6. 

Greek 

f3 Pressure gradient parameter, (4.14). 

7 Fluid specific heat ratio, (3.2). 

7 r Longitudinal intermittancy factor, (4.4). 

6 Displacement thickness, (4.8). 

c Convergence tolerance, (5.2). 

I Ratio of turbulent to molecular viscosity coefficients, (4.13). 

0 Cascade stagger angle, Figure 1. 

9 Momentum thickness, (4.9). 

0 Viscous-layer local to edge temperature ratio, (4.15). 

k Wake streamline curvature, (3.5). 

/z, n T Molecular and effective turbulent viscosities, (4.4). 
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£, T] Cascade axial and “circumferential” Cartesian coordinates, Figure 1; independent 
similarity variables, (4.10). 

p Density, (3.2). 

t w Surface shear stress, e.g. Figure 4(c). 

$ Velocity potential for inviscid flow, (3.1). 

Viscous-layer stream function, (4.11). 

D Flow angle measured from axial direction, positive counter-clockwise, Figure 1. 

oj Relaxation factor, (5.1). 

u7 Total pressure loss parameter, Figure 6. 

Subscripts 

e Viscous-layer edge value, (4.3). 

/, V Inviscid, viscous values, (5.1). 

i Index of solution station, (5.2). 

W, W Wake value, (3.5), (4.22), (4.21). 

TE Trailing-edge value, (4.6), (4.7). 

:poo Far upstream/downstream freestream value of variable, Figure 1. 

Superscripts 

n Global iteration count, (5.1). 

+ , — Upper, lower viscous layer, (4.19), (4.22). 
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Figure 1: Two-dimensional compressor cascade. 
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Figure 2: Semi-Inverse Iteration Procedure. 
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Figure 3: Streamline H-mesh for the EGV cascade. 
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Figure 4: Inviscid ( ) and IVI solutions for EGV cascade at Re = 10 5 ( ) and 

10 6 ( ): (a) pressure coefficient; (b) displacement thickness; (c) surface shear stress. 
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Figure 5: Comparison of IVI ( ) and Navier-Stokes ( ) solutions for the EGV 

cascade at Re = 10 6 : (a) pressure coefficient; (b) displacement thickness; (c) surface shear 
stress. 
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Figure 6: Predictions for the EGV cascade for a range of inlet flow angles: (a) loss parameter; 
(b) exit flow angle; (c) separation point location. 
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Figure 7: Trailing-edge streamline patterns for the EGV cascade: (a) f2_oo = 36 deg; (b) 
fl_oo = 45 deg; (c) = 54 deg. 
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Figure 8: Streamline H-mesh for the HSC cascade. 
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Figure 9: Inviscid ( ) and IVI solutions for the HSC cascade at Re = 10 5 ( ) 

and 10 6 ( ): (a) pressure coefficient; (b) displacement thickness; (c) surface shear 

stress. 
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Figure 11: Streamline H-mesh for the turbine cascade. 
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Figure 12: Results for turbine cascade: (a) comparison of predicted and measured Mach 

number distributions: ( ) IVI; ( ) inviscid; symbols: experimental measurement; 

(b) predicted displacement thickness distribution; (c) predicted surface shear-stress distri- 
bution. 
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Figure 13: Convergence history for the EGV cascade at Re = 10 6 and fi_oo = 40 deg: (a) 
total pressure loss parameter; (b) exit flow angle. 
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A. Details of the Viscous-Layer Solution Procedure 


The viscous-layer equations are solved numerically using an implicit finite-difference ap- 
proach, which is described in this Appendix. The finite-difference approximations used to 
discretize the governing equations and the quasi-linearization applied to the resulting system 
of nonlinear equations is discussed in §A.l. The recursion relations and the associated coef- 
ficients needed to solve the block-tridiagonal system of equations on the surface and along 
the wake are given in §A.2 and §A.3, respectively. 


A.l Finite-Difference Approximations 


The partial derivatives appearing in the governing viscous-layer equations, (4.12) and 
(4.13), are all first derivatives, except for the shear term appearing in the latter equation, 
which introduces a second derivative with respect to //. In the present analysis, first-order 
accurate backward differences are used to approximate d/d£ terms and second-order accurate 
central differences are used to approximate first and second partial derivatives with respect 
to tj. In the description that follows the subscripts i and j are the mesh-point indices in 
the £- and //-directions, respectively, so that the notation (•),- J refers to a quantity evaluated 
at the location (^i,//j). For convenience, the i index will only be used to denote quantities 
evaluated at £ locations other than £,, so that if no i index is given, the variable is evaluated 
at i.e., (-)j refers to (£,, r/j). Surface viscous layers are discretized in the surface normal 

direction using NJ points and 2 NJ — 1 = f NJ 2 points are used across the wake. The index 
j = 1 at blade surfaces (77 = 0 ) and j = NJ at the outer edge of surface viscous layers 
(77 = // e ). In wakes j = 1 at the lower edge (// = //“), j = NJ along the wake reference 
streamline (77 = 0) and j = NJ 2 at the upper edge of the wake viscous layer (77 = 77+). 

Both the direct and inverse modes are used on the blade surface, while only the inverse 
mode is used in the wake; both modes are described below. The procedure used to solve for 
the surface and the wake viscous layers is similar, however, in the latter case, the boundary 
conditions are applied at three different locations, requiring a modification to the surface 
technique, as detailed in §A.3. 

First derivatives of the dependent variables are written as 


and 


where 


and 


fan 

Uv A < 

(A.l) 

3 F\ .. F )+1 - Fj_i 
dy)j~ A7 b+i + Ar h ’ 

(A-2) 

= ti-b- 1, 
A %+ 1 = Vj+1 - V] 


Atjj = T/j - T)j—i . 
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The second derivative with respect to tj is discretized as follows: 


d (o- dF \ ~ 2 

dv \ £ dy ) j ~ 



•^j+i F i 
Aj ij 


— esj-1,2 


Fi-Fj-A 

&Vj - 1 / 


(A.3) 


where, e.g., Cej+1/2 = (^£j+i + p£j)/2. If, as in all of the calculations presented here, 
|1 — Arjj/Arij-i\ <C 1, the ^-derivatives remain formally second-order accurate. 

Before discretizing the equations all variables are written in “delta” form, e.g., Fj = 
Fj + 8Fj , where the superscript g denotes the value of Fj from the previous iteration or 
an initial guess. Thus, rather than solve for the variables themselves, the changes in the 
variables (e.g., 8Fj) are solved for. Following this procedure for the continuity equation, 
(4.12), leads to the expression 


Sfj - 8 fj—i + Pj(8Fj + SFj-i) = Qj , 


(A.4) 


where 


Pj = -Aj/j/2 , 

Qj = —P. ){Pj + Fj- 1) — fj + fj-i ■ 


(A.5) 


In this equation and those that follow, the superscript g has been dropped for convenience. 
The momentum equation can be written, after discretizing and letting a \ = l/(A7/j + Arjj+i) 
and a 2 = 2£,/A£, in the form 


where 


Aj — a i 
Cj = a i 


AjSFj-i + BjSFj + Cj8F j+1 + DjSfj = EijSfi + E 2j , 

2(le) J -i/2 


Aijj 

((&)j+l/2 , (ffij-1/2 


—2a\ I — ' - — I- 

\ Ar/j+i A gj-i 


2(^)i+l/2 


Di 


6iai(l + a 2 )(F j+1 — Fj- 1) , 


+ {/j + a 2 (/j — j)} 


Eu = *i Ff-$j 


and 

Ev 


= —2 aj 


(%l/2n (m ]+ l/2 1/2 \ p (^)i- 1/2 ip 

- 1, Ar,, +1 + A*,- J ^ + A,, ^ 


(A.6) 


— {/? + a 2 (/j “ /i-i.i)} 

- bi [a 2 (2F 3 — + 2 /3Fj] , 


- Fi-u) ~ Mil/; + «2 (fj ~ fi-iMFj+i ~ Fj- i) + (*i*? - *i)P • 


(A.7) 
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The quantity b x controls the FLARE approximation [34] used in separated flow regions; when 
F < 0 at any of the points j — 1, j or j + 1 at either i or i — 1 (except the wall point at 
j = 1), bx is set to zero, otherwise b\ = 1. Note that the term E\S0 is only included in the 
inverse mode, where 0 is one of the unknowns; in the direct mode 80 = 0. 

In deriving the discretized equations, the nonlinear terms, which have the generic form 
UV, are expanded using the “delta” form of the variables as follows: 

UV = ( U 9 + 8U)(V 9 + 8V) = U 9 V a + U 9 8V + V 9 8U + 6U8V . (A.8) 

The governing equations are quasi-linearized by neglecting 0(8 2 ) terms (i.e., 8U8V). The 
truncation error associated with the quasi-linearization is eliminated by performing multiple 
iterations at each solution station. 

The algebraic system of equations resulting from the discretization and quasi-linearization 
of the original partial differential equations is block- tridiagonal in form, allowing standard 
techniques to be used for its solution. The surface and wake solution procedures are simi- 
lar, although the wake analysis is more complicated because of the differences between the 
boundary conditions used for each procedure. 

A. 2 Surface Procedure 

The system of equations applicable to the surface viscous layers is solved using the block- 
tridiagonal inversion procedure described below. The recursion relations for the unknown 
variables 8Fj and Sfj are written as follows: 


SFj = Rj - Sj80 - TjSFj+i 

(A.9) 

Sfj = Lj - MjS0 - NjSF j+1 . 

(A. 10) 


In the direct mode, where 0 is known, the coefficients Sj and Mj are set equal to zero. 

Substituting for SFj - 1 and 6/j_i in the continuity equation (A. 4) using Eqs. (A. 9) and 
(A. 10), solving for Sfj and using this expression and Eq. (A. 9) to eliminate <5Fj_ i yields, 
after collecting terms, the following expressions for Rj , Sj and Tj\ 

Rj = d\[AjRj-\ + Dj(Qj + Z/j-i — PjRj- 1) — Ey] , 

Sj = di[AjSj-! + Dj(Mj_i — PjSj- 1) + Eij] , 

and 

Tj = -d,Cj, (A. 11) 

where 

i , = - B, + DjlNj-, + - 3}.,})] . (A. 12) 
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Returning to the expression for 8 fj obtained by manipulating the continuity equation as 
described in the preceeding paragraph and substituting for 8Fj using Eq. (A. 9) produces the 
following expressions for the coefficients appearing in Eq. (A. 10): 

Lj = Qj + Lj_x — P jRj-i — Rj[Nj- 1 + Pj( 1 — Tj-i)] , 

— PjSj-i — Sj[Nj_i + Pj(l — Tj-i)] 


and 


Nj = -TANj-x + PM ~ Tj-x)} ■ (A. 13) 

The recursion relation coefficients are obtained using Eqs. (A.11)-(A.13), sweeping from 
j = 2 to j = NJ — 1. The values of the coefficients at j = 1 are needed; they are obtained 
from the boundary conditions for / and F, Eqs. (4.17) and (4.18). Since the values of F and 
/ are known at j = 1, we set 6Fj = x = Sfj-x = 0, which gives 

Rx = Sx = Tx = 0 (A. 14) 

and 

L x = Mj = Nx = 0 . (A. 15) 

The solution for 6F : and 6 fj is obtained by sweeping the recursion relations, Eqs. (A. 9) 
and (A. 10), from j = NJ — 1 to j = 2. This requires knowledge of the values of SFpjj, 
6 fwj and 8(3. The quantity SF^j = 0, since Fyvj = 1 (Eq. (4.19)). The technique used to 
determine the value of 6 depends on whether the direct or inverse mode is being used, 
as described below. 

In direct mode the value of /3 is known, hence 8(3 = 0. Writing the continuity equation 
(A. 4) at j = NJ and substituting Eqs. (A. 9) and (A. 10) evaluated at j = NJ — 1 yields 

G/nj = Qnj + Lnj-x ~ PnjRnj-x ■ (A. 16) 


In inverse mode the value of 8fjyj is prescribed using Eq. (4.20), which is rewritten in 
“delta” form as 

772 

tifNJ = h + r) NJ — — y= — /nj . (A. 17) 


Combining this relation, Eq. (A. 4) written at j = NJ and Eqs. (A. 9) and (A. 10) evaluated 
at j = NJ — 1 yields the following solution for 8(3: 


8/3 = 


Qnj — ^Jnj_ + Lnj-x ~ PnjRnj-x 
Mnj-x — PnjSnj-x 


(A. 18) 


With all variables known at j = NJ , Eqs. (A. 9) and (A. 10) are swept from j = NJ — 1 
to j = 2 to solve for 8Fj and 8fj. Finally, the values of Fj and fj are obtained from the 
relations 


Fj = Ff + 8Fj 


(A. 19) 
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and 


fj — ff + &fj 5 (A. 20) 

and the updated value of /? from 

P = P 9 + 6/3 . (A. 21) 

Because the viscous-layer equations are nonlinear, this solution procedure is applied itera- 
tively at each z-station until the solution at £, converges. Convergence occurs when the maxi- 
mum magnitude of the residuals of the continuity and momentum equations, max 2 <j<Arj |Qj| 
and max 2 <j<ATJ-i \F 2 j |> respectively, are less than a specified tolerance, which we usually set 
to 5 x 10 -4 . Convergence is typically obtained in from two to five iterations. 

A. 3 Wake Procedure 

The solution procedure used for the wake viscous layers is very similar to that used for 
the surface layers. However, modifications must be made to account for the application of 

two boundary conditions at the upper and lower edges of the viscous layer (Eq. (4.19)) and 

a third boundary condition along the wake reference streamline (Eq. (4.17)). In addition, 
since the inverse mode is used in the wake, a jump in / from j = 1 (at tj = r}j) to j = N J 2 
(at tj = 77+) is imposed (Eq. (4.22)). 

An additional term is appended to the recursion relations used for the surface viscous- 
layer analysis; the wake recursion relations are given by: 

SFj = Rj - Sj6f3 - TjSFj +1 - UjSfvjt (A.22) 

and 

Sfj = Lj - MjSfl - NjSFj+i - HjSf NJ 2 . (A.23) 

Using the same approach employed for the surface analysis we obtain the coefficients in 
Eq. (A.22). The expressions for /?j, S 3 and Tj are identical to those given in Eq. (A. 11), and 
the value of Uj is given by 


Uj = drlAjUj-1 + D^Hj. r - P 3 U^)\ , (A.24) 

where di is given by Eq. (A. 12). The coefficients Lj, Mj and Nj appearing in Eq. (A.23) are 
defined in Eq. (A. 13), and the additional coefficient Hj is given by 

Hj = Hj-i - PjUj _! - + Pj{\ - Tj-i}) • (A. 25) 

At the lower edge of the viscous layer (tj = T]~ , j = 1) F\ = 1, therefore 8F\ = 0, yielding 

R 1 = S 1 = T 1 =U l = 0 . (A. 26) 

The jump condition for / is rewritten in delta form as 

tifNJ2 ~ ^/i = fl — InJ2 + h + (17+ — V e ) ~ > (A. 27 ) 
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which can be rearranged in the form of Eq. (A. 23) evaluated at j = 1 to give 


Lx = f N J2 - fi - h - (r/ e + - T) e ) + -J= , 

M x = N x = 0 


and 


H\ = — 1 . 


(A. 28) 


With the recursion relation coefficients determined at j = 1, the coefficients at the remaining 
j-locations can be determined. 

The remaining boundary conditions, namely 

6/nj = 0 (A. 29) 

and 

&Fnj2 = 0 , (A. 30) 

must be applied and the value of 8(3 obtained before the system of equations can be inverted 
to obtain 8Fj and 8fj. To accomplish this, the complete system of equations given by the 
recursion relations, the remaining boundary conditions and the continuity equation applied 
at j = N J2 (which has not yet been invoked) are written in matrix form as follows. At each 
point j G [1,7VJ2 — 1] we can write the recursion relations in the form 


where 


and 


AjSU j -|- Bj6U j. j-i + CjSV — dj , 


A? 

Bj 


6U : 

6V 


1 0 

0 1 J ’ 

T, 0 ' 

Nj 0 J ’ 

0 Uj S 3 ‘ 

0 Hj Mj J ’ 
I 

*f: J ’ 

&FnJ 2 
6/nJ 2 

6(3 . 


d i 


R 3 

l 3 


(A.31) 
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The complete matrix system of equations to be inverted has the form 


' Ai 

By 



1 

ei 


6U j 


A2 B 2 



1 

1 

^2 


6 U 2 



Anj 2-2 

B^j 2-2 

1 

1 

eyvj 2-2 


SU_NJ 2-2 




Anj 2-1 

1 

I 

ZNJ2-1 


dU NJ2-1 

. Ci 

C 2 

Cnj 2-2 

Csj 2-1 

1 

9 


6V 


where 


and 


<?; 


Cnj 


Cnj 2-1 


9 


10 0 ' 
0 0 
0 0 

'00* 
0 1 
0 0 


for 1 < j < NJ2 — 2 , j ^ N J , 




PnJ2 — 1 

0 0 

0 0 

•P/VJ2 1 0 

0 0 0 
1 0 0 


d 


QnJ 2 
0 
0 


dr 

d 2 


d_NJ2-2 

dNJ 2-1 

d 

(A. 32) 


Note that e/tfj 2 -i is modified from the standard form of ej so that 


£nj 2-1 = 


Tj Uj Sj 

Nj Hj Mj 


j=NJ 2-1 

Eq. (A. 32) can be rewritten as a partitioned matrix of the form 

r I e \ / 6U 

- I -- 

7 T , 


c 


/ 


6V 


d \ 
K d ) 


(A. 33) 


(A.34) 


The submatrix T is an upper block-bidiagonal matrix — the only nonzero terms are the 
nonzero entries appearing in the 2x2 blocks denoted by A and B in Eq. (A. 32). The row 
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vector C T consists of the 3x1 vector elements C, and represents the continuity equation 
written between NJ 2 and NJ 2 — 1, and the boundary conditions 8 f^j = 0 and 8 F^J 2 = 0, 
respectively. 

Expanding the partitioned matrix and combining the two resulting equations yields 



8 V = (g-C T Y- l e)- x {d-C T Y- l d) . 

(A. 35) 

Letting 

YX = d 

(A. 36) 

and 

YV = e 

(A.37) 

allows Eq. (A. 35) to be rewritten in the form 



8 V = (g-C T Y)-\d-C T X). 

(A. 38) 


Equations (A. 36) and (A. 37) are easily inverted to solve for the elements of X and Y since T 
is an upper block-bidiagonal matrix. Noting that A is the 2x2 identity matrix, the elements 
of X and Y are given by 

X-NJ2-1 — dNJ2-\ i 

Xj = dj - BjXj+i , j = NJ 2 - 2 i 1 


and 


y NJ2-1 


CNJ2-\ i 


Yj = e; - BjY j+l , j = NJ 2 - 2 | 1 . 


(A. 39) 


^r-rr 


The final solution for 6 V is obtained by carrying out the matrix multiplications for C X 
and C Y , forming d — C X and ~g — C Y, inverting the 3x3 matrix g — C Y and solving 


Eq. (A. 38). 

With the elements of 8 V known, the recursion relations, Eqs. (A. 22) and (A. 23), are 
applied to obtain the solution for 8 Fj and 8 fj for 1 < j < NJ 2 — 1. Finally, the values of 
Fj, fj and (3 are obtained using Eqs. (A.19)-(A.21). 
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